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ABSTRACT 

The  growth  of  streamer  trees  in  insulating  fluids  (a  submicrosecond  process  which 
triggers  high-voltage  breakdown)  has  been  simulated  with  a combination  of  parallel- 
coding tools.  Large  grids  and  arrays  display  well  the  multifractal,  self-avoiding 
character  of  the  streamer  trees.  Three  physical  cases  have  been  approximated  by 
different  power-law  weightings  of  the  statistical  growth  filter:  dense  anode  trees, 
in  uniform  field;  sparse  cathode  trees  (a  rarer  experimental  case);  and  ultra-sparse 
anode  trees  (seen  in  some  fluids  of  higher  viscosity). 

The  model  is  contained  in  a software  package  that  is  written  in  Fortran  90 
with  data  parallel  extensions  for  distributed  execution.  These  extensions  encapsu- 
late an  underlying,  invisible  message-passing  environment  (MPI),  thus  enabling  the 
solution  of  memory  intensive  problems  on  a group  of  limited-memory  processors. 
Block-partitioning  creates  processes  of  reasonable  size,  which  operate  in  parallel 
like  small  copies  of  the  original  code.  The  user  needs  only  to  express  his  model  in 
transparent  array-directed  commands;  parallel  interfacing  between  blocks  is  han- 
dled invisibly.  Breakdown  is  performed  in  parallel,  in  each  of  the  local  blocks. 

Results  are  presented  for  experiments  run  on  9 nodes  of  the  IBM  SP2,  and  4 
nodes  of  the  SGI  Onyx,  two  examples  of  multiple-processor  machines.  Display  is 
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carried  out  in  3-D.  Timing  of  the  growth  can  be  shown  by  color  banding  or  by  frame- 
animation  of  the  results.  The  adequacy  of  the  growth  rules  and  size  scaling  are 
tested  by  comparing  the  simulations  against  snapshots  from  high-voltage  discharge 
events. 


2 


1 Introduction 


The  growth  of  fast  streamer  trees  in  liquid-dielectric  insulation  provides 
the  precursor  “leader”  conduction  path  through  which  damaging  flashover 
between  electrodes  can  take  place.  The  high-speed,  variable  nature  of  this 
phenomenon  has  made  its  detailed  mechanism  elusive.  Nonetheless,  a global 
description  of  the  process  may  be  useful  for  its  characterization. 

We  have  applied  stochastic  Laplacian  growth  as  a model  for  feist  streamer 
trees  in  liquid  dielectrics.  Filamentary  dielectric  breakdown  has  been  exten- 
sively discussed  by  Pietronero  and  Wiesmann  [1]  and  others  [2,  3,  4,  5].  Here 
we  construct  a practical  realization  of  the  algorithm  on  a large  Cartesian 
grid,  using  boundary  conditions  which  confine  the  electric  field.  We  exam- 
ine the  effect  of  parameters  (threshold  voltage,  choice  of  power  law)  on  the 
fractal  structure  (dense  or  sparse?)  and  timing  of  the  growth  process. 

The  calculation  of  the  voltage  field  throughout  the  full  volume,  repeated 
after  each  stage  of  breakdown  growth,  is  the  major  computational  burden 
which  calls  for  parallel  methods. 

The  SIMD  (single  instruction  stream,  multiple  data)  model  fits  our  prob- 
lem closely.  We  had  tested  an  earlier  machine- language  version  of  our  algo- 
rithm on  the  CM-2  Connection  Machine  [6].  The  SIMD  version  there  was 
synchronous  and  specific  to  the  machine.  The  present  method  uses  block 
partitioning  under  MPI  [7]  control,  in  a SPMD  (single  program,  multiple 
data)  approach  which  runs  asynchronously.  For  a data-parallel  environment 
which  is  easy  to  use  and  transport,  we  have  written  the  code  in  array-directed 
commands  of  Fortran  90,  together  with  our  library  DPARLIB  [8,  9]  that  en- 
capsulates the  underlying  MPI  [7]  calls.  This  will  run  in  any  environment 
that  is  Fortran  90  and  MPI  enabled,  be  it  networks  of  workstations.  Sym- 
metric Multiprocessors  (SMPs),  parallel  machines,  or  some  combination  of 
these. 

Running  Fortran  90  on  a single  processor  does  not  provide  access  to  the 
memory  that  large  data-parallel  programs  require.  Our  software  has  been 
designed  to  fit  into  the  Fortran  90  syntax  transparently,  but  also  to  run 
across  the  nodes  of  a group  multiprocessor.  It  extends  the  array-directed 
commands  seamlessly,  across  block  boundaries.  For  development  purposes, 
the  number  of  processors  executing  can  be  reduced  to  one,  like  a serial 
version  of  the  code  using  only  plain  Fortran  90  commands.  These  features 
enable  users  to  test  small  versions  of  their  code  on  their  workstations.  Then, 
with  only  minor  changes  in  the  code,  much  larger  versions  can  be  run  on 
machines  like  the  SP-2,  taking  advantage  of  the  larger  memory  and  fast 
communication  network. 
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Physical  Model  and  Algorithm 


The  elements  of  a stochastic  fractal  model  have  been  demonstrated  in  nu- 
merous examples  [10,  11,  12,  13,  14,  15,  16].  The  streamer  tree  is  assumed 
to  be  fully  conductive,  and  electrically  attached  to  the  electrode  from  which 
it  originates.  The  voltage  field  is  defined  over  a set  of  points  in  a rectangular 
grid.  The  streamer  tree  is  a connected  set  of  those  points. 

After  initialization,  the  algorithm  alternates  between  a Laplacian  con- 
vergence procedure  and  statistically  weighted  streamer  tree  growth.  Our 
algorithm  may  be  summarized  as  follows: 


Initialize  voltage  field  and  anode/cathode  configuration 

Set  the  streamer  tree  to  be  the  anode  (including  steirter  needle) 

Run  Gauss-Seidel  over-relaxation  on  Laplace’s  equation  for  200  iterations 

Do  until  streamer  tree  has  reached  cathode 

Set  all  sites  on  streamer  tree  to  voltage  0.0 

Set  all  sites  on  cathode  to  voltage  1.0 

Solve  Laplace’s  equation  via  Gauss-Seidel  over-relaxation 

Find  nearest  neighbors  to  tree 

Do  until  the  set  of  breakdown  sites  is  not  less  than  one 
Set  Breakdown  sites  equal  to  the  empty  set 
Choose,  locally  in  each  block,  red  or  black  (at  random) 

Set  breakdown  candidates  equal  to  neighbours  of  the  color 
Remove  candidates  with  voltage  below  a specified  threshold 
For  each  remaining  site  generate  a weighted  random  number 
Remove  breakdown  sites  whose  voltage  is  below  this  number 

End  Do 

Add  breakdown  sites  to  streamer  tree 
End  Do 


The  two  stages  of  this  algorithm  are  the  solution  of  Laplace’s  equation, 
and  the  growth  of  the  streamer  tree.  We  now  consider  these  in  more  detail. 
Solution  of  Laplace’s  Equation 

Laplace’s  equation  is  fully  solved  throughout  the  interior  region  (filled 
with  dielectric  liquid)  using  the  anode,  starter  needle,  and  streamer  tree 
as  one  boundary  at  zero  potential,  and  the  plane  cathode  as  the  counter- 
electrode. Sides  of  the  cube-volume  support  periodic  boundary  conditions. 
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The  quasi-static  field  approximation  (i.e.,  that  the  voltage  field  always  “catches 
up”  with  the  new  boundary  condition  at  each  stage  of  the  growth)  is  phys- 
ically reasonable  since  breakdown  streamers  are  known  to  advance  at  su- 
personic velocities,  still  slow  by  comparison  with  the  speed  of  light  in  the 
dielectric  liquid. 

A first-order,  six-point-averaging  Laplacian  operator  [17]  on  red-to-black 
and  black-to-red  checkerboard  subgrids  advances  through  Gauss-Seidel  over- 
relaxation. Initially  200  iterations  of  convergence  bring  the  smoothing  preci- 
sion below  one  part  in  ten  thousand,  on  a grid  of  126  x 126  x 126.  (Nominal 
single  iteration  time  is  roughly  10  seconds.)  After  this  first  high-precision 
convergence,  the  epsilon  is  relaxed  (to  .001  - .005)  in  the  interests  of  speed. 
Because  each  growth  stage  creates  only  small  perturbations  on  the  existing 
boundary  conditions,  the  later  cycles  reconverge  rapidly,  in  five  to  ten  itera- 
tions. The  changes  propagate  in  rapidly  diminishing  ripples,  away  from  the 
newest  growth  links  on  the  tree. 

Growth  of  the  Streamer  Tree 

Grid  sites  immediately  adjacent  to  the  tree  are  examined.  If  their  volt- 
ages exceed  a specified  threshold  (or  cutoff)  level,  then  the  surviving  voltages 
are  compared  against  a weighted  distribution  of  random  numbers.  Those 
which  pass  over  the  statistical  hurdle  are  attached  to  the  tree,  and  the  Lapla- 
cian convergence  is  re-cycled  for  the  new,  perturbed  boundary  condition. 
Field  and  growth  stages  alternate  until  the  counter-electrode  is  reached. 

By  contrast  with  some  of  the  stochastic  fractal  models  cited  earlier, 
which  recalculate  Laplace’s  equation  after  each  single  new-growth  site  is 
added,  our  algorithm  considers  simultaneous,  distributed  growth  possibili- 
ties on  all  “red”  or  all  “black”  neighbors  of  the  tree,  during  each  growth 
stage.  Growth  is  favored  near  the  tips  of  the  tree,  where  electrostatic  field 
lines  converge  strongly,  and  the  voltage  gradients  are  largest.  While  the 
details  of  the  mechanism  for  fast  streamer  growths  are  still  to  be  deter- 
mined from  experiment,  we  ask,  “To  what  laws  of  shape  and  timing  does 
the  observed  growth  conform?” 

Weighting  of  the  growth-probability  power  law,  as  a function  of  electric 
field  strength,  determines  the  degree  of  bushiness  (or  fractal  dimension)  of 
the  final  structure.  At  one  extreme,  Sanchez  et  al.  [18]  have  shown  that  a 
fourth  power  (or  higher)  dependence  on  field  strength  produces  a single,  self- 
avoiding strand  which  resembles  a directed  random  walk.  Side-branching  is 
minimized.  At  the  other  limit,  Witten  and  Sander  [19]  showed  that  a linear 
filtering  against  uniformly-distributed  random  numbers  results  in  a dense 
growth  form,  which  is  also  found  in  diffusion-limited  aggregation,  from  a 
dififerent  mechanism.  In  this  linear  regime,  a mass  of  twigs  advances  with  a 
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“front”  of  thickly-spaced  growing  tips.  A recent  review  by  Erzan,  Pietronero, 
and  Vespignani  [20]  considers  ways  of  characterizing  fractal  dimension  in  this 
class  of  problems. 

In  Figure  1 we  see  a log-log  plot  of  breakdown  probability  (for  a sin- 
gle link  between  grid  vertices)  versus  power  of  the  electric  field  strength, 
measured  as  the  voltage  difference  between  the  conducting  tree,  at  zero  po- 
tential, and  its  immediate  neighbors  on  the  Cartesian  grid.  The  power  law 
determines  the  likelihood  of  growth  from  the  tip,  as  against  side  branching. 
Recall  that  Laplace’s  equation  is  size-scalable  without  limit.  By  adopting  a 
power-law  probability,  we  may  expect  self-similar  growth  over  large  ranges  of 
both  size  and  field  strength.  This  is  not  viewed  as  a microscopic  description, 
but  rather  as  a law  controlling  the  global  growth  form. 

Notice  how  miniscule  the  breakdown  probability  becomes  for  cube-law 
and  higher  powers,  on  the  left  half  of  the  graph.  The  range  shown,  of  field 
strengths,  is  that  encountered  in  growth  across  a 126  x 126  x 126  grid  vol- 
ume. As  the  streamer  tip  advances  from  anode  to  cathode,  the  diminishing 
gap  leads  to  progressively  higher  fields  around  the  leading  branches,  while 
electrostatic  screening  causes  the  fields  around  the  base  of  the  tree  to  drop 
back  below  threshold,  so  that  growth  is  cut  oflF  there.  Because  the  electric 
field  is  confined  by  limiting  boundary  conditions,  and  the  growth  process 
starts  and  stops  at  an  abrupt  threshold,  we  do  not  expect  a “pure”  fractal 
dimension  throughout  the  growth.  The  tree  is  sparser  as  it  leaves  the  con- 
centrating needle  on  the  primary  electrode,  and  denser  as  it  approaches  the 
counter-electrode. 

The  algorithm  must  operate  uniformly  across  a wide  range  of  growth 
rates.  For  square-law  and  higher  powers,  the  starting  probability  near  the 
anode  needle  (where  the  field  strength  only  slightly  exceeds  the  threshold 
value  of  0.015  to  0.10)  is  very  small  - tens  or  hundreds  of  empty  statisti- 
cal tries  may  occur  between  single  breakdown  events.  By  contrast,  linear 
weighting  leads  to  tens  of  events  per  growth  stage,  in  the  shorter  gaps,  ap- 
proaching the  cathode. 

The  uniform  treatment  of  slow  and  fast  stochastic  processes  assumes 
that  the  discretized  physical  model  approximates  a Markov-process  master 
equation  [21,  22],  whose  time  variable  is  proportional  to  the  number  of  trial 
instances. 
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3 Implementation 


Our  code  was  first  implemented  in  a serial  version,  on  Sun  and  SGI  work- 
stations, and  was  run  interactively  on  small  grids.  This  enabled  extensive 
testing  with  visual  tracking  of  detailed  printouts  to  the  screen.  This  serial 
code  was  implemented  in  standard  Fortran  90  and  exploits  that  language’s 
array-oriented  syntax  and  intrinsic  functions.  As  examples  of  the  powerful 
simplifications  available,  we  describe  some  of  the  details. 

The  code  makes  extensive  use  of  array  arithmetic,  the  WHERE  con- 
struct, array  functions  such  as  CSHIFT,  MAXVAL,  COUNT,  and  the  dy- 
namic memory  allocation  features  of  Fortran  90.  Our  program  logic  is  ex- 
pressed through  combinations  of  real  and  logical  arrays,  the  WHERE  as- 
signment mask  of  Fortan  90,  and  a selection  of  spreading  and  collecting 
operations. 

For  example,  our  Laplacian-operator  subroutine  makes  vigorous  use  of 
Fortran  90’s  CSHIFT  operation.  The  sum  of  neighbors  calculation  was 
implemented  in  the  serial  version  with  the  following  code: 


biksum  = cshift 

(blkphi,  +1, 

1) 

+ cshift 

(blkphi , 

-1. 

1) 

& 

+ cshift 

(blkphi,  +1, 

2) 

+ cshift 

(blkphi , 

-1, 

2) 

& 

+ cshift 

(blkphi,  +1, 

3) 

+ cshift 

(blkphi , 

-1, 

3) 

CSHIFT  is  the  standard  Fortran  90  intrinsic  function  that  performs  cir- 
cular shifts.  The  array  (biksum)  that  results  from  this  calculation  con- 
tains at  each  grid  point  the  sum  of  the  six  neighbors  of  that  grid  point  in 
the  original  array,  using  periodic  boundary  conditions.  The  division  of  the 
three-dimensional  grid  into  alternating/adjacent  red  and  black  subgrids  is 
a feature  brought  over  from  statistical  mechanics  simulations.  Black-to-red, 
and  red-to-black  averaging  steps  are  carried  out  across  all  processes,  at  each 
loop.  This  routine  is  global  over  the  entire  internal  voltage-field  array. 

We  also  made  extensive  use  of  the  “modules”  feature  of  Fortran  90.  This 
simplified  the  debugging  process  by  enabling  the  compiler  to  catch  certain 
types  of  coding  errors.  In  particular,  the  compiler  checked  that  all  subroutine 
calls  had  properly  constructed  calling  sequences. 

The  statistical  selection  was  coded  by  the  following  steps.  By  contrast 
with  the  Laplacian  routine,  these  are  local  operations,  executed  indepen- 
dently in  each  block. 

1.  Find  all  nearest  neighbor  sites  to  the  growth  tree,  within  the  process 
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block.  This  procedure  produces  a logical  mask  array  that  is  .true,  at 
the  neighbor  sites  and  .false,  elsewhere. 

2.  Select,  at  random  for  each  block,  either  the  red  or  black  neighbors 
within  that  block.  The  mask  array  of  neighbors  is  modified  to  re- 
flect this  choice.  Since  we  are  considering  simultaneous  but  indepen- 
dent growth  across  the  full  distributed  array  of  neighbor  sites,  this 
red/black  choice  discourages  selection  of  adjacent  sites,  which  might 
lead  to  growth  concentration  or  cause  a locally  enhanced  perturbation 
of  the  voltage  field. 

3.  Allocate  a one-dimensional  REAL  array  with  one  element  for  each 
surviving  neighbor.  If  a power  law  higher  than  linear  is  desired,  then 
additional  arrays  of  the  same  length  are  allocated. 

4.  Fill  the  one-dimensional  array(s)  with  uniformly  distributed  random 
numbers.  The  DPARLIB  dp_uni  is  used  to  generate  random  numbers 
that  are  independent,  across  all  of  the  parallel  processes. 

5.  Calculate  the  MAX  array  of  these  random  arrays.  The  resulting  array 
contains  a weighted  distribution  at  one  power  lower  than  the  desired 
breakdown  power  law.  Step  7 below,  the  statistical  filtering,  is  equiv- 
alent to  a single  order  of  integration  for  the  distribution  of  voltages 
which  will  survive  this  hurdle. 

6.  UNPACK  the  MAX  array  through  the  neighbor  mask  array.  This 
creates  a three-dimensional  array  with  the  weighted  random  numbers 
at  the  surviving  neighbor  sites. 

7.  Compare  the  voltages  to  the  weighted  random  array  at  the  neighbor 
sites.  WHERE  they  pass,  add  new  sites  to  the  existing  growth  tree. 
Write  out  their  grid  positions,  trial  numbers,  and  field  strengths. 

8.  If  no  sites  have  been  added  to  the  tree,  return  to  Step  2. 

9.  If  a site  has  been  added  to  the  tree,  return  to  the  Laplacian-convergence 
stage. 

Note  that  the  approach  used  in  steps  3 to  6 substantially  reduces  memory 
and  computational  requirements  by  generating  the  weighted  distribution  of 
random  numbers  in  minimally-sized  arrays.  The  resulting  linear  array  is 
UNPACKed  into  the  full  size  three-dimensional  array  only  when  necessary. 
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The  Fortran  90  array  instructions  which  we  needed  were:  WHERE, 
CSHIFT,  UNPACK,  MAXVAL,  MAX,  MOD,  COUNT,  ANY,  ALLOCATE, 
DEALLOCATE.  In  the  block-parallel  mode,  some  are  used  in  local  fashion 
by  the  individual  processor,  addressing  only  the  portions  of  the  major  arrays 
seen  within  its  own  block. 

4 Conversion  to  Parallel  Operation 

The  DPARLIB  operations,  which  were  required  as  substitutions  or  additions 
to  enable  parallel  operation,  were  dp_cshift,  dp_uni,  dp_count,  dp_any, 
dp_synch  [8,  9].  In  addition,  the  dpJnitialize  and  dp_query Jayout 
routines  were  called  at  the  start  of  the  run,  to  set  up  the  division  of  domains 
between  processors. 

The  adaptation  of  the  serial  code  was  accomplished  through  the  use  of 
the  DPARLIB  subroutine  library.  This  is  a Fortran  90  library  developed  at 
NIST  to  facilitate  the  transition  of  serial  applications  to  a parallel  computing 
environment  and  the  development  of  new  parallel  programs.  DPARLIB  is 
built  on  the  Message  Passing  Interface  (MPI),  an  industry-standard  library 
for  passing  data  and  coordinating  the  activities  of  multiple  processes  in  a 
parallel  operating  environment. 

DPARLIB  is  designed  to  be  used  in  a single-program-multiple-data  (SPMD) 
programming  approach.  In  other  words,  multiple  copies  of  the  same  pro- 
gram are  running  simultaneously,  and  each  copy  is  processing  a different 
portion  of  the  data.  In  particular,  DPARLIB  provides  simple  mechanisms 
to  divide  very  large  arrays  into  blocks,  each  of  which  is  handled  by  a sep- 
arate copy  of  the  program.  In  practice,  this  means  that  the  researcher  can 
write  parallel  code  that  looks  almost  identical  to  serial  code.  In  our  case, 
the  code  could  be  written  as  though  addressed  to  a single  active  grid-node, 
and  its  immediate  neighbors.  Fortran  90,  extended  across  block  boundaries 
by  DPARLIB,  executed  each  instruction  on  all  2 million  sites  of  each  array. 

In  this  way,  an  existing  serial  program  can  often  be  converted  to  a parallel 
program  with  few  changes  by  using  DPARLIB,  which  plays  the  role  of  a 
high-level  language  for  block  parallelism. 

DPARLIB ’s  emphasis  on  array  handling  is  designed  to  mesh  with  Fortran 
90’s  array  syntax  and  intrinsic  array-handling  functions.  Much  of  DPARLIB 
consists  of  parallel  versions  of  the  intrinsic  array  functions  such  as  CSHIFT 
and  MAXVAL.  Because  DPARLIB  is  coded  entirely  in  standard  Fortran  90 
and  depends  only  on  MPI,  it  is  portable  to  any  environment  that  provides 
those  two  resources. 


9 


The  adaptation  of  the  dielectric  breakdown  code  was  straightforward 
and  involved  the  following  steps: 

- Add  calls  to  DPARLIB  housekeeping  functions  for  initialization  and 
specification  of  the  mapping  of  data  arrays  to  an  array  of  processes.  (dpJnitialize 
and  dp.query -layout) 

- Make  all  distributed  data  arrays  allocatable,  and  allocate  the  arrays 
based  on  information  provided  by  DPARLIB  calls. 

- Change  various  Fortran  90  intrinsic  calls  to  the  corresponding  DPAR- 
LIB calls.  This  simply  involves  changing  the  name  of  the  called  routine; 
e.g.,  cshift  becomes  dp_cshift  and  maxval  becomes  dp_m£ixval. 

Once  this  conversion  was  complete,  each  process  ran  the  same  program 
with  its  own  block  of  data.  It  used  no  explicit  reference  to  data  within  other 
blocks,  except  as  information  was  received  through  collecting,  spreading, 
and  shifting  operations  in  DPARLIB.  The  parallel  version  could  be  run  at 
once,  on  a network  of  workstations  with  the  LAM  environment  of  MPI  [23]. 

Our  Laplacian  operator  is  first-order,  and  requires  only  the  exchange  of 
data  for  single  layers,  between  domain  blocks  [25]. 

A refreshing  feature  of  DPARLIB  is  that  communication  between  pro- 
cesses is  handled  automatically.  Thus,  the  typical  boundary-layer  transfers 
of  data,  as  described  by  Gropp,  Lusk,  and  Skjellum  [7]  do  not  require  ex- 
plicit consideration  by  the  code  writer.  dp_cshift  includes  - invisibly  - the 
necessary  “shadow”  sites  and  cross-boundary  transfers  of  data  which  must 
take  place  at  the  planes  between  process  domains  [25,  17].  The  user  can 
ignore  this  complex  and  error-prone  aspect  of  code  writing. 

There  were,  however,  two  crucial  points  in  this  program  at  which  we 
had  to  take  some  care  about  the  parallel  execution  of  multiple  copies  of  the 
code:  writing  the  streamer  tree  sites  to  a disk  file,  and  testing  for  whether 
the  streamer  tree  had  reached  the  cathode.  At  these  points,  the  application 
code  must  explicitly  deal  with  the  fact  that  multiple  copies  are  executing  in 
parallel. 

At  the  end  of  the  tree  growth  part  of  the  algorithm,  the  grid  sites  that 
have  been  added  to  the  streamer  tree  are  written  to  a data  file.  Because  these 
sites  will  be  found  in  multiple  processes  we  have  these  multiple  processes 
trying  to  write  to  the  same  output  file  at  the  same  time.  This  will  result 
in  errors  on  many  systems,  including  the  SUN,  SGI,  and  IBM  systems  on 
which  we  were  running  these  tests. 

This  problem  can  be  solved  by  having  each  process  write  to  its  own  out- 
put file.  These  separate  output  files  could  then  be  merged  after  the  program 
is  complete.  Where  multiprocessor  access  to  a common  file  is  possible,  we 
instead  implemented  a “round-robin”  solution  that  involves  the  synchroniza- 
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tion  of  the  processes  using  the  DPARLIB  routine  dp_sync.  In  this  solution, 
each  process  in  turn  opens  the  common  output  file,  appends  the  necessary 
information,  then  closes  the  file,  so  that  no  two  processes  are  writing  to  the 
output  file  at  the  same  time. 

The  other  point  in  the  algorithm  at  which  the  code  must  explicitly  deal 
with  the  parallel  execution  of  the  code  is  in  testing  whether  the  streamer 
tree  has  reached  the  cathode.  This  is  because  only  a portion  of  the  tree,  and 
of  the  cathode  may  be  present  in  each  of  the  processes.  The  issue  is  handled 
by  a simple  test  based  on  the  which  part  of  the  full  data  array  has  been 
assigned  to  each  process.  If  a process  contains  a portion  of  the  cathode, 
then  the  appropriate  part  of  the  streamer  tree  array  that  resides  in  that 
process  participates  in  a test  of  whether  the  tree  has  reached  the  cathode. 

5 Results  and  Interpretation 

The  crucial  test  was  “Can  streamer  growth  proceed  seamlessly  through 
the  defined  block  boundaries?”  Figure  2 illustrates  an  example  of  a linear- 
weighted  streamer  simulation,  grown  with  nine  process  blocks.  The  block 
partition  boundaries  were  crossed  as  expected,  by  the  DPARLIB  version 
of  the  program  which  contained  no  explicit  reference  to  boundaries,  in  the 
array-directed  commands.  Neither  the  Laplacian  nor  the  statistical-test  sub- 
routine faltered  in  crossing  - meaning  that  data  transfer  had  proceeded 
successfully. 

Each  run  delivers  an  output  list  of  vertices,  recorded  with  the  index 
number  of  the  statistical  (red-black)  trial.  Knowing  that  they  form  a singly- 
connected  tree,  we  find  the  links  between  sites  by  means  of  a backward- 
search  algorithm.  Because  the  statistical  trial  number  - the  clock-tick  for 
“Monte  Carlo  time”,  in  a stochastic  process  [24]  - has  been  recorded  for 
each  breakdown  event,  a simulation  of  time  progression  in  the  global  growth 
is  possible. 

Display  is  carried  out  in  three  dimensions,  using  color  banding  to  mark 
the  time  history  of  growth.  Animation  of  frames  is  also  readily  achieved. 
Visual  presentation  of  the  data  can  be  compared  against  high-speed  pho- 
tographs from  experiments  [26,  2],  which  are  two-dimensional  shadowgraphs 
taken  in  side  view.  (In  examining  the  results  of  very  large,  dense  fractal 
growths,  the  combination  of  color  banding  and  three-dimensional  dynamic 
rotation  can  often  call  attention  to  features  of  the  growth  which  are  not 
obvious  from  static  projections.) 

We  include  two  examples  of  photography  from  experiment  [26].  Figure 
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3 (a  and  b)  shows  two  separate  samples  of  dense  anode  streamer  formation, 
each  growing  in  less  than  30  nanoseconds  in  n-hexane,  in  a “uniform”  field 
configuration  (meaning  that  the  growth  commences  close  to  the  anode,  in 
a small  gap  between  large-radius  spheres).  The  time  sequence  proceeds 
from  top  to  bottom,  at  100-nanosecond  intervals.  The  third  frame  in  each 
sequence  captures  the  start  of  the  “leader”  breakdown,  with  the  central 
shock  wave  growing  in  diameter  in  the  last  two  frames. 

Figure  2 is  one  Cartesian  projection  of  a three-dimensional  dense  growth, 
whose  isometric  view  is  displayed  in  Figure  4.  Figure  5 shows  the  orthog- 
onal 2-D  projection.  What  we  note  with  some  interest  is  the  well-defined 
cone-shaped  upper  envelope  to  the  growth;  Stricklett  et  al.  [26]  have  called 
particular  attention  to  this  feature  of  the  experimental  observations.  This 
cone  represents  the  boundary  at  which  further  lateral  growth  is  cut  off  by 
screening,  causing  the  electric  field  to  drop  below  the  cutoff  (threshold) 
value  for  discharge.  This  feature  of  the  numerical  experiment  is  undergoing 
further  investigation. 

For  this  case  of  growth  with  a uniform  (linear)  statistical  filter,  the  ma- 
jor divisions  of  the  tree  have  branched  laterally  to  fill  the  intervening  spaces 
with  twigs.  They  remain  self- avoiding,  however.  The  heavily  weighted  track 
in  Figure  4 (and  later  in  Figure  8)  is  the  first  continuous  connecting  path  be- 
tween electrodes  - it  simulates  the  ionization  path  along  which  the  “leader” 
flashover  will  occur. 

Figure  6 displays  a sparse  fast  cathode  streamer,  formed  under  similar 
conditions,  but  in  the  reverse  polarity.  The  positioning  of  the  electrodes  has 
been  reversed;  the  cathode  is  at  the  top,  the  anode  at  the  bottom.  Here 
the  timing  sequence  of  frames  proceeds  firom  bottom  to  top.  In  this  Figure 
6 case,  the  fast  “secondary”  cathode  streamer  begins  from  a primary  slow- 
growing  bushy  structure;  Watson  [27]  and  Fenimore  [28]  have  described  the 
primary  stage  by  a bubble  simulation.  (The  bipolarity  of  the  anode  and 
cathode  fast  streamer  processes  presents  a further  challenge  for  physical 
interpretation  - how  is  filamentary  tree  growth  possible  in  both  directions 
under  like  conditions?) 

Figures  7 and  8 show  growth  with  a square-law  statistical  filter,  for 
comparison  against  the  fast  cathode  streamer  in  Figure  6.  In  this  case,  we 
note  that  both  experiment  and  simulation  show  a single  major  “trunk” , or 
column  of  growth,  proceeding  half  way  across  the  gap,  then  dividing  into  3 
or  4 major  branches,  whose  envelope  lies  within  a more  narrow  angle  than 
the  dense  anode  streamers  of  Figure  3a  and  3b.  Here  the  forward- directed 
field  is  shaping  the  growth.  The  lateral  twig  growth  from  the  major  branches 
has  not  filled  the  intervening  gaps;  screening  does  not  dominate. 
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Figure  9 provides  an  example  of  successive  growth  frames,  as  a function 
of  cumulative  statistical  trial  numbers.  Development  of  the  growth  against 
“Monte  Carlo  time”  shows  notable  acceleration  as  the  gap  is  shortened  and 
the  resultant  electric  field  (voltage  gradient)  increases.  In  fact,  the  lower 
half  of  the  growth  takes  place  in  roughly  20  percent  of  the  total  elapsed 
duration. 

The  fronts,  or  growth-tip  envelopes  for  both  power-law  weightings  show 
considerable  variability  from  run  to  run.  This  is  not  surprising,  if  we  recall 
that  the  statistical  filters  are  weighted  in  favor  of  extreme  (high)  values  of 
electric  field,  producing  a runaway  tendency.  Note  that  the  “leader”  con- 
duction path  from  needle  to  counter-electrode,  shown  in  heavy  line  weight, 
resembles  a directed  random  walk.  In  some  instances,  two  branches  of  the 
leader  may  reach  the  counter-electrode  on  the  same  statistical  try. 

6 Time-compression  for  higher  powers 

The  extremely  low  probabilities  shown  in  the  lower  half  of  Figure  1 compli- 
cate the  trial  of  cube-law  and  higher  powers.  For  low  values  of  initial  field 
strength,  many  thousands  of  empty  statistical  trials  would  be  required  to 
secure  a very  few  steps  of  growth.  Our  solution  is  to  re-weight  the  statis- 
tical test,  after  the  manner  of  the  BKL  algorithm  [24],  so  that  the  a-priori 
likelihood  of  at  least  one  discharge-event  is  increased  to  one. 

The  changed  routine  for  higher  powers  proceeds  as  follows,  after  step  3 
in  the  sequence  which  was  described  earlier: 

• PACK  the  surviving  neighbor- voltage  values  onto  the  first  linear  real 
array. 

• square,  cube,  or  otherwise  exponentiate  this  array  onto  a second  linear 
array.  Sum  the  values  of  the  elements  in  this  array,  to  a constant 
we  call  ‘sumsq’.  It  represents  the  a-priori  probability  of  at  least  one 
discharge-event. 

• Fill  a third  linear  array,  of  the  same  length,  with  uniformly  random 
numbers.  Multiply  these  by  the  coeflBcient  ‘sumsq’  - this  is  equivalent 
to  “compressing”  Monte  Carlo  time. 

• WHERE  the  second  array  fails  to  exceed  the  random  array,  weighted 
by  the  coefficient  ‘sumsq’,  set  the  values  of  the  first  linear  array  equal 
to  zero 
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• UNPACK  the  revised  first  array  through  the  position  mask  of  the 
original  neighbors 

• WHERE  the  unpacked  3-D  array  is  now  greater  than  zero,  add  new 
sites  to  the  existing  tree. 

• Update  the  statistical  trial  number  (Monte  Carlo  time  tick)  by  an 
integer  equal  to  FLOOR(1.0/sumsq).  Thus  we  retain  an  estimate, 
rather  than  a count,  of  the  number  of  empty  trials  which  would  secure 
one  successful  discharge. 

• Return  to  the  Laplacian-convergence  subroutine. 

Again,  Fortran  90  has  furnished  powerful  array-directed  operations  for 
composing  the  program  logic.  The  additional  computing  burden  is  small, 
since  these  operations  are  performed  on  few  elements  - numbering  only  in 
the  thousands. 

The  modified  algorithm  has  been  successfully  run  for  cube-law  and  fourth- 
power  examples.  Figure  11  illustrates  simulation  of  a cube- law  streamer, 
run  on  a 128  x 128  x 128  grid,  with  a short  anode  needle,  2 grid  steps  in 
length.  The  spareness,  restricted  lateral  branching,  and  directed  alignment 
are  characteristic  of  growth  with  this  exponent  in  the  power  law. 

Recent  reports  by  Miyano  and  collaborators  [29]  display  fast  anode 
streamers  in  D40  perfluoro  polyether  and  FCIO  perfluoro  carbon.  Their 
photographs  show  strong  similarities  to  Figure  11. 

Thus,  we  note  that  physical  counterparts  to  our  power-law  simulations 
with  exponents  1,  2 and  3 are  found  in  the  experimental  observations  on 
fast  filamentary  streamers. 

7 Limitations 

There  are  several  clear  limits  in  the  present  realization  of  the  model: 

• Diagonal  breakdown  paths  between  grid  nodes  have  not  been  included 
in  the  simulation.  Increased  flexibility  of  directional  choice,  at  each 
event,  would  be  a step  towards  greater  physical  realism,  albeit  impos- 
ing a greater  bookkeeping  burden. 

• The  Laplacian-convergence  accuracy  has  been  relaxed  in  a manner 
which  may  tend  to  favor  tip-growth  versus  sidebranching,  to  a modest 
extent.  This  was  considered  acceptable  at  the  present  semi-quantitative 
condition  of  comparison  against  experiment. 
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• Visual  comparison  against  experiment  is  affected  by  details  of  the  vi- 
sual perception,  rendering  of  line-width  weight,  and  3-D  obscuration. 
(Thus,  comparing  Figure  2 with  Figure  4,  we  see  that  the  projection 
along  the  Cartesian  axis  looks  distinctly  more  transparent.) 


8 Summary 

The  completed  program  is  a robust  and  flexible  tool  for  investigations  of 
physics.  All  communications  and  coordination  between  copies  of  the  pro- 
gram are  hidden  within  calls  to  DPARLIB  routines.  The  code  is  expressed 
and  compiled  in  standard  Fortran  90,  using  DPARLIB  routines  that  closely 
emulate  standard  Fortran  90  intrinsic  functions. 

The  MPI  complexity  is  invisible.  Thus,  DPARLIB  has  fllled  in  the 
needed  elements,  for  a high-level  language  of  block  parallelism. 

The  use  of  parallel  methods  has  aflForded  improvements  in  physical  re- 
alism: large  array  domains,  parallel  breakdown  across  the  entire  tree  struc- 
ture, and  clocking  of  Monte  Carlo  time.  Many  problems  in  the  physical 
sciences  are  well-adapted  to  such  a SPMD  parallel  treatment.  Large  arrays 
corresponding  to  a spatial  domain  can  be  partitioned  across  many  processes. 
Such  algorithms  can  often  be  simply  expressed  in  Fortran  90,  and  DPARLIB 
enables  a quick  transition  from  a serial  to  a parallel  environment,  without 
the  need  to  learn  complex  communications  techniques  between  processes. 

For  portability,  the  code  was  run  under  SGI  (MIPSpro),  and  IBM  (xpf) 
compilers.  The  size  of  the  executable  code  depends  on  the  choice  of  grid 
bounds.  The  executable  version  on  the  IBM  SP-2,  for  a grid  of  126  x 126  x 
126,  was  roughly  15  megabytes  for  each  of  9 processes.  Speed  was  not 
the  major  objective;  rather,  the  ability  to  handle  a large  number  of  very- 
large-sized  arrays  within  the  local  process  memory.  Often  a salient  issue 
in  physics  computing  is  space,  not  time.  To  easily  spread  a problem  over 
many  compute  nodes  can  mean  the  difference  between  being  able  to  study 
a problem,  and  not  being  able  to  do  it  at  all. 

Comparative  timings  have  been  carried  out  for  the  routine  running  on 
the  IBM  SP2  and  on  two  SGI  multiprocessors  (Onyx  and  Origin).  These 
are  for  the  linear  statistical  comparison,  running  on  a (small)  64  x 64  x 64 
lattice.  Because  the  choice  of  random  number  seeds  varies  with  the  number 
of  processes,  the  actual  executions  are  for  slightly  different  problems.  The 
numbers  give  a general  idea  of  speed  on  dedicated  multiprocessors. 

The  timings,  in  seconds,  using  1,  2,  4,  and  8 processes  are  as  follows: 

• IBM  SP2:  12445,  5336,  3256,  1770 
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• SGI  Power  ONYX:  7767,  3590,  1588,  824 

• SGI  Origin:  5467,  2310,  1328,  690 

Full  run  times  on  the  larger  128  x 128  x 128  lattice  are  roughly  ten 
times  as  long,  since  the  larger  problem  involves  a longer  path  from  anode 
to  cathode,  as  well  as  greatly  enlarged  Laplacian  convergence  tasks.  The 
execution  of  Figure  11,  for  example,  was  carried  out  on  four  processors  of  a 
shared  SGI  Power  Onyx,  with  overnight  turnaround  time. 

DPARLIB  is  based  only  on  standard  Fortran  90  [30]  and  on  the  Message 
Passing  Interface  (MPI)[7].  MPI  is  an  industry  standard  library  that  imple- 
ments message  passing  in  parallel  computing  environments.  Both  Fortran 
90  and  MPI  are  available  on  a variety  of  computers  and  work  stations.  LAM 
is  an  open  cluster  enviroment  for  work  stations  or  multiprocessors  [23]. 

Our  dielectric-breakdown  code,  named  CADMUS,  will  be  available  through 
our  web  site 

http : //www. itl .nist .gov/div895/ sasg/ 
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10  Disclaimer 


Certain  commercial  equipment  and  software  may  be  identified  in  order  to 
adequately  specify  or  describe  the  subject  matter  of  this  work.  In  no  case 
does  such  identification  imply  recommendation  or  endorsement  by  the  Na- 
tional Institute  of  Standards  and  Technology,  nor  does  it  imply  that  the 
equipment  or  software  is  necessarily  the  best  available  for  the  purpose. 


17 


References 

[1]  L.  Pietronero  and  H.J.  Wiesmann,  Z.  Phys.  B70,  87  (1988) 

[2]  W.G.  Chadband,  IEEE  Trans.  Elect.  Ins.  v.  23,  no.4  697  (1988). 

[3]  W.G.  Chadband  and  T.M.  Sufian,  IEEE  Trans.  Elect.  Ins.  v.20,  no. 2, 
239  (1985). 

[4]  S.  Satpathy,  in  Fractals  in  Physics^  Elsevier  Science  Publishers  B.V.,  p. 
173  (1986). 

[5]  T.  J.  Lewis,  in  The  Liquid  State  and  its  Electrical  Properties^  E.E.  Kun- 
hardt,L.G.  Christophorou,  and  L.E.  Luessen  eds  (Proc.  NATO  Adv. 
Study  Inst.  July  5-17,  1989,  Sintra,  Portugal),  N.Y.  (Plenum)  1988, 
pp.  431-453. 

[6]  We  are  indebted  to  Larry  Davis  and  Mitchell  Murphy  of  the  Univer- 
sity of  Maryland  Institute  for  Advanced  Computing  Systems,  for  their 
assistance  and  for  access  to  the  CM-2. 

[7]  W.  Gropp,  E.  Lusk,  and  A.  Skjellum,  Using  MPI,  Cambridge  (MIT 
Press),  1994,  esp  pp.  60-80. 

[8]  Judith  E.  Devaney  and  John  Hagedorn,  “Transforming  a MIMD  Hard- 
ware Environment  into  a SIMD  Programming  Environment” , National 
Institute  of  Standards  and  Technology  internal  paper.  May  7,  1996. 

[9]  John  Hagedorn  and  Alan  Heckert,  “DPARLIB  User’s  Guide”,  draft  of 
NIST  software  library  documentation,  March  18,  1997. 

[10]  L.  Niemyer,  L.  Pietronero,  and  H.J.  Wiesmann,  Phys.  Rev.  Lett.  52, 
1002  (1984). 

[11]  L.  Pietronero  and  H.J.  Wiesmann,  J.  Stat.  Phys.  36,  909(1984). 

[12]  H.J.  Wiesmann  and  L.  Pietronero,  in  Fractals  in  Physics,  Elsevier  Sci- 
ence Publishers  B.V.,  p. 151(1986). 

[13]  L.  Pietronero,  C.  Evertsz,  and  H.J.  Wiesmann,  ibid.,  p. 159(1986). 

[14]  M.  Murat,  ibid.,  p.  169(1986). 

[15]  E.  Arian,  P.  Alstrom,  A.  Aharony,  and  H.E.  Stanley,  Phys.  Rev.  Lett. 
63,  2005(1989). 


18 


[16]  P.D.  Beale  and  P.M.  Duxbury,  Phys.  Rev.  B37,  2785  (1988). 

[17]  Jagdish  J.  Modi,  Parallel  Algorithms  and  Matrix  Computation,  Claren- 
don Press,  Oxford  (1988),  pp. 135-143. 

[18]  A.  Sanchez,  F.  Guinea,  L.M.  Sander,  V.  Hakim,  and  E.  Louis,  Phys. 
Rev.  E48,  no. 2, 1296  (1993).  In  these  treatments,  idealizations  such  as 
radial  growth  to  a remote  spherical-shell  boundary,  with  normalized 
probabilities,  favor  growth  having  a simple,  well-defined  fractal  dimen- 
sion. 

[19]  T.A.  Witten  and  L.M.  Sander,  Phys.  Rev.  Lett.  47,  1400  (1981);  Phys 
Rev.  B27,  5686  (1983). 

[20]  A.  Erzan,  L.  Pietronero,  and  A.  Vespignani,  Rev.  Mod.  Phys.  67,  no. 3, 
545-604  (1995). 

[21]  M.M.  Clark,  L.M.  Ruff,  and  H.L.  Scott,  Computers  in  Physics  10,  no. 6, 
584-590  (1996). 

[22]  K.  Binder  and  D.W.  Heermann,  Monte  Carlo  Simulation  in  Statistical 
Physics.  Springer,  Berlin,  1988,  esp  pp.  28-29. 

[23]  D.  Burns  and  R.B.  Daoud,  “Lam;  an  open  cluster  environment  for  mpi” , 
in  Supercomputing  Symposium  ’94,  June  1994,  Toronto,  Canada.  Code 
available  at  http//www. osc.edu/lam.html. 

[24]  I.Beichl  and  F.E.  Sullivan,  “(Monte  Carlo)  Time  after  Time”,  Com- 
putational Science  and  Engineering  4,  no. 3 (July-Sept.  1997), pp. 91-94. 
See  also  A.B.  Bortz,  M.H.  Kalos,  and  J.L.  Lebowitz,  J.  Computational 
Physics  17,10  (1975)  esp.  Appendix  B,  “Computation  of  Stochastic 
Time  Variable”,  also  J.G.  Amar  and  F.  Family,  Phys.  Rev.  Lett.  74, 
2006  (1995). 

[25]  For  a similar  treatment  of  a second-order  operator,  see  A.  Villareal  and 
J.A.  Scale,  Computers  in  Physics  11,  no.4,  388-399  (1997),  esp.  p.  392. 

[26]  K.L.  Stricklett,  E.F.  Kelley,  H.  Yamashita,  H.  Kawai,  and  C.  Fenimore 
in  Conference  Record,  IEEE  Int’l  Symposium  on  Electrical  Insulation, 
1990,  p.  161-164.  These  data  were  taken  by  means  of  a high-speed  “op- 
tical storage  line”  technique  which  limited  the  available  spatial  resolu- 
tion. We  are  greatly  indebted  to  Dr.  Stricklett  for  making  his  original 
data  available  to  us. 


19 


[27]  P.K.  Watson,  in  Conference  Record  of  the  10th  International  Confer- 
ence on  Conduction  and  Breakdown  in  Dielectric  Liquids,  Grenoble, 
France,  Sept.  10-14,  1990,  pp.  22-28. 

[28]  C.  Fenimore,  in  Conf.  Record  of  the  1988  IEEE  Int.  Symposium  on 
Electrical  Insulation,  Cambridge,  MA,  June  5-8,  1988,  pp. 27-30. 

[29]  K.  Miyano,  K.  Yamazawa,  M.  Pompili,  C.  Mazzetti,  and  H.  Yamashita, 
in  1997  Annual  Report  of  the  Conference  on  Electrical  Insulation  and 
Dielectric  Phenomena,  Minneapolis,  Oct. 19-22,  1997,  v.  II,  pp.  640- 
643  (available  from  IEEE  Service  Center,  Piscataway,  N.J.).  We  refer 
especially  to  Figures  3 and  4 in  this  paper.  The  authors  associate  the 
properties  of  the  observed  streamers  with  the  viscosity  of  the  perfluoro 
liquids  used. 

[30]  J.C.  Adams,  W.S.  Brainerd,  J.T.  Martin,  B.T.  Smith,  and  J.L.  Wa- 
gener,  Fortran  90  Handbook^  McGraw-Hill,  New  York,  1992. 


20 


11  Figure  Captions 


Figure  1.  Log- log  plot  of  breakdown  probability  (discharge  probability  for 
an  individual  link)  versus  electric  field  strength,  for  power  laws  with  integral 
exponents.  The  horizontal  scale  is  electric  field  strength  at  the  candidate 
links,  in  the  discretized  step  approximation  (difference  between  the  zero 
voltage  of  the  streamer  tree  and  the  voltages  at  the  neighbor  grid  nodes). 
The  limit  value  1.0,  to  which  all  fields  are  normalized,  is  the  value  which 
would  be  reached  between  flat-plate  electrodes  separated  by  one  grid  step. 

Figure  2.  Lateral  projection  of  a linear-weighted  streamer  simulation 
on  a grid  of  126  x 126  x 126,  showing  the  boundary  planes  between  the 
nine  process  blocks.  The  centered  starting-needle  extends  15  grid  units 
from  the  upper  (anode)  plane;  cathode  plane  is  at  the  bottom.  Periodic 
boundary  conditions  are  applied  at  the  side  walls  of  the  cube.  In  this  run, 
597  statistical  tries  produced  27203  discharged  grid  links.  Cutoff  (threshold) 
voltage  was  set  at  0.05.  Real  time  for  the  calculation  on  the  IBM  SP-2  was 
5:17:58. 

Figure  3 (a  and  b).  two  examples  of  fast  anode  streamer  discharges, 
grown  in  n- hexane  in  less  than  50  nanoseconds  [26].  The  gap  between  elec- 
trodes is  3 mm;  radius  of  the  spherical  electrodes  is  1.27  cm.  Anode  at 
top,  cathode  at  bottom.  Peak  voltage  was  250  kV.  Time  sequence  is  top 
to  bottom,  in  100-nanosecond  intervals.  The  growth  occurs  in  less  than 
30  nanoseconds,  before  the  third  frame  shown  in  each  sequence.  The  dark 
volume  is  filled  with  fine  filaments,  incompletely  resolved  here.  Other  ex- 
periments have  shown  the  filaments  to  have  individual  diameters  as  small 
as  4 micrometers.  They  are  composed  of  individual  small  straight  linear 
segments,  most  being  of  length  30  micrometers  or  longer. 

Figure  4.  Isometric  projection  of  the  3-D  dense  anode-streamer  simu- 
lation in  Figure  2.  The  branching  structure  in  3-D,  which  was  hidden  in 
the  Figure  2 projection,  becomes  evident.  The  tree  is  dense,  but  the  ma- 
jor branches  remain  self-avoiding  because  the  electric  field  is  screened  out 
of  the  gaps  between.  Scaled  to  the  dimensions  of  Figure  4,  the  individual 
discharged  links  (each  one  grid  interval)  would  have  length  approximately 
30  micrometers.  Apparent  visual  density  is  dependent  on  the  choice  of  line 
width  in  the  graphical  reconstruction.  Compare  against  Figure  3.  The 
leader  path  connecting  the  electrodes  is  shown  in  heavy  line  weight. 

Figure  5.  2-D  projection  of  the  simulation  in  Figure  4,  orthogonal  to 
Figure  2.  Note  the  well-defined  cone  angle  of  the  upper  envelope,  a feature 
corresponding  to  experiment. 
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Figure  6.  Fast  (secondary)  cathode  streamer,  produced  under  condi- 
tions similar  to  Figure  3.  Cathode  at  top;  anode  at  bottom.  The  time 
sequence  now  proceeds  from  bottom  to  top.  The  tree  structure  is  notice- 
ably sparser.  This  event  occurs  with  distinctly  lower  probability  than  the 
fast  anode  streamer.  The  bushy  growth  attached  to  the  cathode  at  top  is 
the  slower  (primary)  streamer.  Both  anode  and  cathode  fast  streamers  pro- 
vide the  “leader”  breakdown  path  which  leads  to  high  current  flow  and  the 
formation  of  the  shock  wave. 

Figure  7.  Lateral  projection  of  a square- law- weighted  streamer  simu- 
lation on  a grid  of  126  x 126  x 126.  4227  statistical  tries  produced  4432 
discharged  links.  Threshold  voltage  was  set  at  0.015.  To  compensate  for  the 
very  low  starting  probability,  all  voltage  values  on  neighbor  sites  were  raised 
by  0.020  when  comparing  against  the  statistical  filter.  Real  time  13:32:37. 

Figure  8.  Isometric  view  of  Figure  7.  Compare  against  Figure  6.  The 
cathode  is  now  at  the  top;  anode  at  the  bottom. 

Figure  9.  Color-banded  version  of  the  streamer  tree  shown  in  Figure  8. 
The  color  palette  is  a linear  representation  of  cumulative  statistical  tries,  a 
measure  of  Monte  Carlo  time. 

Figure  10.  Animated  screen  development  of  Figure  8,  against  cumulative 
statistical  tries  (Monte  Carlo  time).  The  acceleration  of  growth  across  the 
second  half  of  the  gap  is  notable.  Note  that  the  duration  of  the  full  growth 
would  correspond  to  a period  of  around  30  nanoseconds,  or  less  than  the 
interval  between  frames  in  Figure  3 and  Figure  6. 

Figmre  11.  Cube-law  streamer  simulation  on  a grid  of  128  x 128  x 128. 
The  concentrating  needle  is  two  grid  steps  in  length.  Spare  appearance  of 
the  filament  tree,  limited  lateral  branching,  and  directed  alignment  along 
the  uniform  field  are  characteristic  of  growth  simulations  with  this  power- 
law  exponent.  This  simulation  was  produced  with  the  time-compression 
version  of  the  algorithm.  It  resembles  experimental  observations  by  Miyano 
et  al.  [29] 
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Figure  1:  Log-log  plot  of  breakdown  probability  versus  electric  field  strength. 
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Figure  2;  Lateral  projection  of  a linear-weighted  streamer  simulation  on  a 
grid  of  126  X 126  X 126,  showing  the  boundary  planes  between  the  nine 
process  blocks. 
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Figure  3:  Two  examples  of  fast  anode  streamer  discharges.  Anode  at  top, 
cathode  at  bottom. 
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Figure  4;  Isometric  projection  of  the  3-D  structure  in  Fig.  2.  The  leader 
path  is  shown  in  heavy  line  weight. 
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Figure  5:  2-D  projection  of  Figure  4,  orthogonal  to  Figure  2. 
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Figure  6:  Fast  (secondary)  cathode  streamer,  produced  under  conditions 
similar  to  Figure  3.  Cathode  at  top;  andode  at  bottom. 
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Figure  7:  Lateral  projection  of  a square-weighted  streamer  simulation  on 
grid  of  126  X 126  X 126. 
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Figure  8:  Isometric  view  of  Figure  7.  Compare  against  Figure  6.  The 
cathode  is  now  at  the  top;  anode  at  the  bottom. 
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Figure  9:  (Original  showed  color-banding  of  the  streamer  tree  in  Figure  8, 
to  represent  statistical  tries  as  a measure  of  Monte  Carlo  time.) 
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Figure  10:  Animated  screen  development  of  Figure  8,  against  cumulative 
statistical  tries  (Monte  Carlo  time). 
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Figure  11:  Cube-law  streamer  simulation  on  a grid  of  128  X 128  X 128. 
Spare  appearance  of  the  filament  tree,  limited  lateral  branching,  and  directed 
alignment  along  the  uniform  field  are  characteristic  of  growth  simulations 
with  this  power-law  exponent. 
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